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Abstract 



Aims. Stars twinkle because their light propagates through the atmosphere. The same phenomenon is expected on a longer time scale 
when the light of remote stars crosses an interstellar turbulent molecular cloud, but it has never been observed at optical wavelengths. 
The aim of the study described in this paper is to fully simulate the scintillation process, starting from the molecular cloud description 
as a fractal object, ending with the simulations of fluctuating stellar light curves. 

Methods. Fast Fourier transforms are first used to simulate fractal clouds. Then, the illumination pattern resulting from the crossing 
of background star light through these refractive clouds is calculated from a Fresnel integral that also uses fast Fourier transform 
techniques. Regularisation procedure and computing limitations are discussed, along with the effect of spatial and temporal coherency 
(source size and wavelength passband). 

Results. We quantify the expected modulation index of stellar light curves as a function of the turbulence strength -characterised by 
the diffraction radius Rdiff- and the projected source size, introduce the timing aspects, and establish connections between the light 
curve observables and the refractive cloud. We extend our discussion to clouds with different structure functions from Kolmogorov- 
type turbulence. 

Conclusions. Our study confirms that current telescopes of ~ 4 m with fast-readout, wide-field detectors have the capability of 
discovering the first interstellar optical scintillation effects. We also show that this effect should be unambiguously distinguished 
from any other type of variability through the observation of desynchronised light curves, simultaneously measured by two distant 
telescopes. 
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1. Introduction 



Galaxy: disk - Galaxy: halo - Galaxy: structure - Galaxy: local interstellar matter - ISM: 



This paper is a companion paper to the observational results 
published in Habi bi et al.l (201 lb), and it focusses on the sim- 
ulation of the scintillation effects that were searched for. Cold 
transparent molecular clouds are one of the last possible can- 
didates for the missing baryons of cosmic structures on differ- 



ent scales (jPfenniger & Combes 1994 Pfenniger & Revaz 2005 



and [McGaugh et al. 2010[ >. Since these hypothesised clouds do 



not emit or absorb light, they are invisible for the terrestrial ob- 
server, so we have to investigate indirect detection techniques. 
Our proposal for detecting such transparent clouds is to search 
for the scintillation of the stars located behind the transparent 
medium, caused by the turbulence of the cloud (IMoniez 20031 
and IHabibi et al.l 201 lb). The objective of this technical paper 
is to describe the way we can connect observations to scintil- 
lation parameters through a realistic simulation. We used these 
connections in the companion paper (Habi bTet al.1 201 lb) to es- 
tablish constraints both from null results (towards SMC) and 
from observations pointing to a possible scintillation effect (to- 
wards nebula B68). Similar studies of propagation through a 
stochastic medium followed by Fresnel diffraction have been 



made by Coles et al. (1995) and for use in radio-astronomy by 
|Hamidouche and Lestrade (2007)| 

We first introduce the notations and the formalism in sec- 
tion [2] Then we describe the different stages of the simulation 
pipeline up to the production of simulated light curves in section 



[3] We study the observables that can be extracted from the light 
curve of a scintillating star, and in particular, we check the ex- 
pected modulation amplitude properties in section|4] The discus- 
sion is extended to non-Kolmogorv turbulence cases in section 
|5] In section|6]we use the results from the simulation pipeline to 
optimise the observational strategy for discovering scintillating 
stars, and indicate some perspectives in the conclusion. 

Complementary information on observations made with 
the ESO-NTT telescope and on the analysis based on the 
present simulations are to be found in our companion paper 
(IHabibi et al.l 2011b). 



2. Basic definitions and formalism 

The formalism described in this section has been inspired and 
adapted from the radioastronomy studies ( jNarayanj 1992). But 
at optical wavelength, the scintillation is primarily due to the 
refraction through dense clouds of H2 + He instead of the in- 
teraction with the ionised interstellar medium. The origin of the 
stochastic phase fluctuation experienced by the electromagnetic 
wave when crossing the refractive medium is the phase excess 
induced by the stochastic fluctuation of the column density due 
to the turbulence ([Moniez 20Q3J: 



(2tt) 2 a 

(f>(x u yi) = — - — Nl(xi,yi) 



(1) 
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Figure 2. The phase-delay variations near the average for a 
simulated refractive screen with N x x N y — 20, 000 x 20, 000 
pixels, Ai = 22.6 km, and Rdiff =100 km. The grey scale ranges 
between +50 X 27T rad ( clear regions correspond to an excess of 
phase with respect to the average). The zoom (inset) illustrates 
the self- similarity of the simulated screen (grey scale amplitude 
of 5 x 2n rad). 
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Figure 3. Phase structure functions D^r) for a phase screen 
with Rdiff = 500 km. Blue line is the initial (theoretical) struc- 
ture function. Red line is reconstructed from one of the realisa- 
tions of the phase screen through simulation. The black curve 
is obtained from the numerical integration of the initial phase 
spectral density sampled as in the simulation. 



Notations 




Figure 1. Geometric configuration. The source is located in 
the (x2,yi) plane, the screen contains the refractive structure, 
and the observer is located in the (xo,yo) plane. A\(xi,y\) and 
A\ (x\ , yi) are the amplitudes before and after screen crossing. 



where x\ and y\ are the coordinates in the cloud's plane, per- 
pendicular to the sightline (see figure [TJ. Here $(x\,y\) is the 
phase delay induced to the wavefront after crossing the cloud, 
Nl(x\,yi) is the cloud column density of Ho molecules plus He 
atoms along the line of sight, a is the medium polarisability, and 
A the wavelength. The phase delay here scales with A~ l ,'m con- 
trast to the radioastronomy where it scales with A. Since other 
sources of phase delay, such as the geometrical delay induced by 
scattering from cloud inhomogeneities, are negligible, the thin 
screen approximation can be used, and the cloud can be con- 
sidered as a 2D scattering screen whose optical properties are 
mapped by the phase screen <p(xi,yi). The statistical properties 
of the phase screen are described by the phase structure function 



D<t>(xi,yi). By assuming an isotropic turbulence ( |Narayanj l992): 
D^xuyi) = D$(r) 



= ([^(xi+x'^yi+yl)- (j)(x[,y[)\ )( x [, y [) = 



Rdi 



'If 



1-2 



(2) 



where the first expression is averaged over the plane positions 
(Xj, jj), r = yjx^ + y 2 v B is the turbulence exponent -equals 1 1/3 

for Kolmogorov turbulence- and the diffraction radius, Rdiff, 
is the transverse distance on the phase screen where the phase 
changes by one radian on average. The diffraction radius can be 
expressed in terms of the cloud parameters (Habibi et al. 201 lb); 
assuming Kolmogorov turbulence it is given by 



Rdi ff (A) = 263km 



lpm 



IQAU 



IQAU 



0~3„ 



10 9 c 



,-3 



(3) 



where L z is the cloud size along the sightline, L ou , is the tur- 
bulence outer scale, and 1x3,, is the dispersion of the volumic 
number density in the mediurrQ. Here we assume the gas to be 
a mixing of H2/He with 24% He by mass (corresponding to the 
primordial abundances) and therefore < a >= 0.720 x 10~ 30 m 3 . 
In this expression, the cloud parameters are scaled to the val- 
ues given by the Pfenniger-Combes model for the clumpuscules 
(the tiniest cloudlets of the molecular cloud). In the NIR band, 
the diffraction radius of a typical clumpuscule is expected to be 
Rdiff ~ 500 km. 



1 A cloud column of width L z can include several turbulent struc- 
tures with outer scale L m „. The direct relation of L m „ with the turbulence 
strengh explains why /?#// increases in conjunction with this parame- 
ter. By contrast, since the column density increases with L-, then the 
refraction also increases, thus decreasing Rdiff. 
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The phase statistics of the screen can be equivalently de- 
scribed in Fourier space by the phase spectral density: 



S4,(q x ,q y ) = S 



d2 

R diff 



2(27^-1/08) 



{Rdiff q)' 



(4) 



where Fourier coordinates q x and q y have inverse length dimen- 
sion, q = ^jq 2 x + q\, and f(j3) = 2 ^(^j^ 2 ' i s a constant. 

After crossing the cloud, the distorted wavefront of a point 
source propagates toward the observer and produces an illumi- 
nation pattern on the observer's plane given by 



Io(x Q ,y ) = -f/i(x ,yo), 



(5) 



where Io(xo,yo) is the light intensity on the observer's plane, L s 
is the source luminosity, z\ is the source-screen distance (see 
Fig.[TJ, and h(xo,yo) is given by the Fresnel-Huygens diffraction 
integral after considering the Fresnel and the stationary phase 
approximations (IBorn & Wolf 20021 IMoniez 20031 : 



h(x ,yo) = 



2nR\ £[ 



^ (A-Q--T 1 )-+0 Q -V 1 )- 



dx\dy\ 



(6) 



where Rf = y/Aza/27T is the Fresnel radius^ and zo is the screen- 
observer distance. The Fresnel radius can be expressed as 



2214 km 



lpm 



Co 



\kpc 



(7) 



At the typical distance of a halo object (~ 10 kpc), Rf ~ 
7000 km for NIR wavelengths. Because of the motion of the 
cloud with respect to the Earth-source line-of-sight, the illumi- 
nation pattern sweeps the observer plane, so that a terrestrial ob- 
server receives fluctuating intensity light from the point source. 
This effect, the scintillation, has two different scattering regimes 
dUscinski 19771 Tatarskii & Zavorotnyi 1980 >, a weak regime 



(Rdiff > Rf) and a strong regime (Rdiff < Rf)- In the present 
studies, we concentrate on the strong regime, which clearly is 
easier to detect, but some realistic configurations may involve 



the intermediate regime studied by Goodman & Narayan (2006) 
For the strong regime, there are two different modes of flux 
variations (|Narayan 1992, see also lRickett 19861 Rumsey 1975 



ISieber 1982b . The first one is the diffractive mode with length 
scale corresponding to the screen's scale of phase variations 
Rdiff given by equation (01. The resulting "speckles", with typ- 
ical size of the order of Rdiff, are shown in figure |4] The corre- 
sponding time scale of the light fluctuations is tdiff = Rdiff /Vt' 



tdiffW 



(8) 



2.6s 



lpm 



10AU 



10AU 



10\ 



,-3 



100 km /s 



where Vj is the sightline relative transverse motion. Therefore 
fast flux variations are expected with a typical time scale of 
tdiff ~ few s. The second variation mode is the refractive mode 
associated to the longer length scale called refraction radius: 



RrefW = 



tea 



RdiffW 



30,860/tm 



lpm 



1 kpc 



RdiffW 



1000 km 



■ (9) 



2 This definition assumes that zo « Zi- In the general case, z l 
should be replaced by Zq 1 + z\ l . 



This natural length scale corresponds to the size, in the ob- 
server's plane, of the diffraction spot from a patch of Rdiff{A) in 
the screen's plane. This is also the size of the region in the screen 
where most of the scattered light seen at a given observer's po- 
sition originates. Our convention for R re f differs from |Narayan| 
(1992) by a factor 2n since it emerges naturally from the Fourier 
transform we use for calculating the illumination pattern (see 
formula[T4l. and it also matches the long distance-scale flux vari- 
ations visible in figure |4] The corresponding time scale is given 

by t re f = R re flV T : 



t r ef(A) — 5.2 min. 



lpm 



1 kpc 



Rdiff{X) 



1000km 



Vt 



100 km/ s 



-i 



(10) 



3. Simulation description 



In this section, we describe the simulation pipeline with some 
numerical tricks, from the generation of the phase screen in- 
duced by turbulent gas, up to the light versus time curves ex- 
pected from realistic stars seen through this gas. The steps in 
this pipeline are listed below: 



- The simulation of the refractive medium: in subsection 13.11 
we describe the generation of a phase screen and examine the 
impact of the limitations caused by the sampling and by the 
finite size of the screen by comparing the initial (theoretical) 
and reconstructed diffraction radii. 

- The computation of the illumination pattern: in subsection 
13.21 we first describe the calculation of the illumination pat- 
tern produced on Earth by a point monochromatic source as 
seen through the refractive medium. We explain the tech- 
nique for avoiding numerical (diffraction) artefacts caused 
by the borders of the simulated screen, and discuss the crite- 
rion on the maximum pixel size to avoid aliasing effects. The 
pattern computation is then generalised to extended poly- 
chromatic sources. 

- The light curve simulation: we describe in subsection !3.3l the 
simulation of the light fluctuations with time at a given po- 
sition induced by the motion of the refractive medium with 
respect to the line of sight. 

3. 1. Simulation of the phase screen 

Numerical realisations of the 2D phase screen — made of N x xN y 
squared pixels of size Ai — are randomly generated from the 
phase spectral density S,f,(q x ,qy), determined by the choice of 
Rdiff in relation (0). Such a phase screen — with the desired sta- 
tistical properties — is obtained from the random realisation of 
a Fourier transform F^,(q x , q y ) in a way that makes the ensemble 
of such realisations satisfy the relation 



< \Ft/>(qx, q y 



> realisations' 



L x L y S^(q= yjqj+qj), (11) 



where L x and L v are the screen physical size, and the average 
is the ensemble averaging over the different realisations. For 
each (q' x ,q J y ) vector associated to a pixel (i,j), we generate a 
random complex number F^{q' x ,q{) = (f^ + i.fjfy/ V2 where 
each component f£ and f£ J spans the Gaussian distributiorQ of 

zero mean, and L x L y S ^(q = ^jq' x 2 + q{ 2 ) width. In this way, re- 
lation QT| is automatically satisfied when averaging on a large 



The turbulence is considered as a Gaussian field 
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Figure 4. Typical illumination pattern from a point-like source. Figure 5. Typical illumination pattern from an extended source 

Here Rdiff =100 km, the screen is at zo = 160 pc, A = 2.16pm, produced through the same screen as in Fig. @ with source ra- 

then Rf= 1300 km and R re f = 106000 km. The typical length dius r s = 0.5/? Q , located at zo + Z\ = Ikpc + 160pc (R s 

scale of the small-scale speckles is Rdiff, and the scale of the 53, 000 km). The small-scale speckles are smeared, and only 

larger structures is R re f- The white square shows our fiducial the larger scale fluctuations survive. The white square shows 

zone with a margin of L, n = R re f/2 from the borders. Grey the restricted fiducial zone with margin of L,„ = R re f/2 + R s 

scale range from to 4 times the mean intensity. The image from the borders. Grey scale ranges +20% around the mean 

has 20, 000 x 20, 000 pixels, each with a 22.6 km side. intensity. 



number of such realisations, since the average of \F$(q x , q y )\ 2 
for the ensemble of (q x ,q Y ) vectors with the same module q 
equals L x L y S,p(q) by construction. The phase screen <p{x\,y{) is 
finally obtained by numerically computing the inverse Fourier 
transform of this phase spectrum F,p{q x , q v ). Figure [2] shows a 
phase screen generated by assuming Kolmogorov turbulence (J3 
= 11/3). 

3.1.1. Preliminary checks, limitations 

To check the accuracy of the numerically generated phase screen 
(Fig. we recomputed the phase structure function D™ c (r) 
(and consequently R^te) from the generated phase Fourier trans- 
form F$(q x , q y ), and we compared it with the theoretical phase 
structure function (Eq. First, the spectral density is recom- 
puted from the generated F^,(q x , q y ) using relation 



ST(q) = 



i\F4q x ,q y )\ 2 ) 



+</; 



L X Ly 



(12) 



where the average is performed on the (q x , q y ) coordinate^ span- 
ning the circle of radius q = ^jq 2 + q 2 . The corresponding phase 
auto-correlation function is then given by Fourier transform: 



(q)e M ^ r dq, 



4 Here again, we infer the ergodicity property that allows us to re- 
place the ensemble averaging with an average on directions from only 
one realisation. 



% rec (r) = J*™ JT qS l <p ec {q)e lniqrcas0 d8dq 
= f 2nq S ™ c {q) J (27rqr) dq, 

Jdrni* 



(13) 



where Jq is the Bessel function. The recomputed structure func- 
tion is then given by D' ec (r) = 2(% rec (0) - % rec (r)), and the value 
of R™ is deduced from its definition D™ c (R™° ff ) = 1. 

In figure [3] we show the theoretical phase structure func- 
tion of a turbulent medium with Rdiff = 500 km -for which 
D<p(r = 500 km) = 1 by definition-, and the recomputed (effec- 
tive) structure function from one of the realisations of the screen. 
From this recomputed function, we find R^f ~ 540 km since 
D', ec {r x 540 km) = 1 . To find the origin of the difference with 
the input value Rdiff = 500 km, we replaced S + ec {q) in equation 
( fT3l by the theoretical spectrum S$(q) sampled as in the simu- 
lation (number of pixels N x x N y ~ 14, 000 x 14, 000 with pixel 
size Aj = 28.85 km). Then we computed the integral ( fT3l l numer- 
ically with the same integration limits {q m i„, q max )^ as the simu- 
lation. The integration result differs only by a few percent from 
the function recomputed from the simulated screen. We showed 
that the black curve approaches the blue curve when q m i„ — > 
and q max — > oo. This means that the sampling is mainly responsi- 
ble for the difference between and D£ ec . Since our simulation 
is limited by the number of pixels, we lose the contributions of 
the large and small scales in the recomputed Rdiff - The only way 
to push back this limitation is to generate larger screens (larger 



In ID: q m 



and q m 



2A[ • 
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N x and N y ) with higher resolutions (smaller Ai) to cover wider 
interval of spatial frequencies which in return needs higher com- 
putational capacities (see also Sect. l3.4t . 

3.2. Illumination pattern 



L 











L - 2L m 






-< >- 






Fiducial Zone 











Figure 6. The smoothing function S F(x). L is the screen size, L, 
is the margin from the screen borders. 



where L m = \0Rf is the margin length from the borders of 
the screen with size L. We multiply the function G{x\,y\) by 
SF(xi) x SF(y\) in expression ( TBI . Since the Fresnel integral 
is dominated by the contribution of the integrand within a disk 
that is within a few Fresnel radii, it is sufficient to smooth the 
discontinuity of G{x\,y\) within a distance of a few Fresnel 
radii (here L„, = lORf). We tested the efficiency of this reg- 
ularisation procedure by checking that the simulated illumina- 
tion pattern from a point-source projected through a uniform 
phase screen was -as expected- also uniform beyond our pre- 
cision requirements (< 1%) (Habibi 20 1 lab . To define a reliable 
fiducial domain in the observer plane excluding the regions that 
are only partially illuminated owing to the diaphragm, we also 
delimited an R re //2 margin from the borders of the illumina- 
tion pattern, corresponding to the typical radius (half-size) of the 
large-scale luminous spots (see also |Coles et al. (1995)| l. Figure 
|4] shows the pattern produced by a point source through a tur- 
bulent medium with Rdiff = 100 km located at z\ = 160 pc 
from the Earth at wavelength A = 2.162/im. Since the corre- 
sponding Fresnel radius Rf — 1300 km is larger than the diffrac- 
tion radius, the regime is the strong scintillation regime. The hot 
speckles (of typical size Rdiff ~ 100 few) can be distinguished 
from the larger dark/luminous structures that have a typical size 
of R ref ~ 2nR 2 F IR diff * 100, 000 Am 



To obtain the illumination pattern on the observer plane, we 
numerically compute the integral © which can be written as a 
Fourier transform: 



h(x ,y ) = 



where 



1 



2nRj 



G(x\,y\) = exp 



\FT [G(xi,vO] 



0(xi,yi) + 



v * Mil J> InRl 



2R 2 F 



2 \ 



(14) 



(15) 



Here, (xuyi) and (xo,yt>) are the screen and observer coordi- 
nates, respectively. Coordinates (f x ,f y ) are the conjugated vari- 
ables in Fourier space. Before computing expression (IT4l . a reg- 
ularisation procedure for G{x\,y\) has been defined to avoid 
computational artefacts. 

3.2.1. Screen regularisation 

Since integral (TBI is computed numerically, the coordinates 
(xi,y{) have discrete (integer) values describing pixel position 
centres on the screen, to allow simple combinations of illumi- 
nation patterns with different pixel sizes (corresponding to dif- 
ferent wavelengths). That the integration domain is limited is 
physically equivalent to computing the Fresnel integral within 
a diaphragm with the size of the screen. In this case, we face 
a parasitic effect: the light diffraction from the sharp edges of 
the diaphragm. This causes rapid intensity variations at the bor- 
ders of the observer plane. To attenuate this effect and remove 
the resulting diffraction fringes, we multiply the screen intensity 
transmission by a 2D smoothing function. We define the follow- 
ing ID smoothing function SF(x) (see Fig- [6j: 

' i [i + sin( s _ |)j o < x < L m , 

1 J—itji X -Zj l—iyyt . 



SF(x) = 



i [l + sin( 7r( ^ L '" ) + f )] L - L m < x < L, 



otherwise, 



3.2.2. Effect of sampling 

Here, we discuss some limitations caused by the pixellisation. 
The screen should be sampled often enough to avoid the aliasing 
effects. Aliasing happens when G(xi,y\) contains frequencies 
that are higher than the Nyquist frequency fy yq = 1 / (2Ai), where 
Ai is the pixel size. In relation ( tlSt . G contains two length scales, 
the diffraction and the Fresnel radii {Rdiff and Rf). Rdiff is the 
characteristic length of the phase screen variations 0(xi,yi); it is 
at least necessary that A; < Rdiff 12 in order to sample phase 
variations up to 1 / 'Rdiff spatial frequency. Rf appears in the 

2Rj 

erates as x\ and yi increase, and aliasing occurs as soon as the 
distance between two consecutive peaks is smaller than 2A\. In 
one dimension we therefore expect aliasing if 



quadratic term exp 



the oscillation of this term accel- 



(jci +2A 1 ) 2 -j^ 



2R\ 



i.e. — > n 
Ai 



Rf_ 
Ai 



> 2tt, 



- 1, 



(16) 



where X\/A\ is the distance to the optical axis, expressed in pix- 
els. The condition Ai = Rf/2 would be obviously insufficient 
here to avoid aliasing, since beyond ~ 1 1 pixels only from the 
optical axis (x\/Ai > 11 pixel), the quadratic term would be 
under-sampled. In practice, for the configurations considered in 
this paper, the Fresnel radius is > 1000 fez. When Ai = 22km 
in the simulation, the aliasing starts at xi/A[ > 6490 pixels 
from the centre of the image, corresponding to more than 140/? f, 
which is large enough to cover the sensitive domain related to the 
stationary phase approximation. 

3.2.3. Extended source (spatial coherency) 

The illumination pattern of a scintillating extended source is 
given by the convolution product of the illumination pattern 
of the point-like source with the projected source limb profile 
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dMoniez 2003l & IHabibi 20"TTib . 



hxt = —P, * K 



(17) 



where L s is the source luminosity, and the normalised limb pro- 
file is described as a uniform disk: 



P r (xo,yo) 







x^+yl<R s 



otherwise, 



where R s — & r s is the projected source radius on the observer's 
plane, and r s is the source radius. Figure [5]shows the convolution 
of the pattern of figure |4] with the projected profile of a star with 
radius r s = Q.5R Q , located at zo (160pc) +z\ {\kpc)=\ .16 kpc 
(R s = 53,000 km). High-frequency fluctuations due to diffrac- 
tive speckle disappear, and the pattern loses contrast, but the 
variations on R re f scale remain visible. As the convolution in- 
volves a disk of radius R s , we can perform the calculation only 
at a distance greater than R s from the borders. We therefore de- 
fine a new fiducial zone by excluding a margin of R re //2 + R s 
from the initial borders. Any statistical analysis will be made 
within this zone to be safe from any border perturbation. 

3.2.4. Polychromatic source (time coherency) 

The illumination patterns shown in Figures |4]and[5]are computed 
for a monochromatic source (fixed A), but observations are done 
through filters with finite-width passbands. To take the contri- 
butions of different wavelengths to the pattern into account, we 
superimpose the illumination patterns obtained with the same re- 
fractive structure (the same column density fluctuations) at dif- 
ferent wavelengths 0. We have considered the passband of the 
SOFI camera in K s band and approximated it as a rectangular 
function over the transmitted wavelengths with a central value 
2. 162 fim and width 0.275 fim. Twenty-one illumination patterns 
were computed for 21 regularly spaced wavelengths within the 
[2.08 , 2.28] yum interval, and were co-added to simulate the il- 
lumination pattern through the K s passband. We checked that 
the spacing between successive wavelengths was small enough 
to produce a co-added image with a realistic residual modula- 
tion index, by studying this index as a function of the number 
of monochromatic components equally spaced within the full 
bandwidth. We found that an asymptotic value is reached with 
a co-added image made up of only about ten components. 

Figure [7] (left and centre) shows a comparison between 
monochromatic (up) and K s passband (down) illumination pat- 
terns of a point-like source. The speckle pattern is attenuated 
when the light is not monochromatic (or equivalently when time 
coherency is limited). This is caused by the small decorrelation 
bandwidth 5Ad ec of the strong diffractive scintillation regime, 
which is given by 

6A dec /A~(R diff (A)/R F f (18) 

( Narayan| 1992, lGwinn et al~9 98). This chromatic effect results 
from the high sensitivity of the constructive interference con- 
dition with the wavelengtfQ. In the case of Fig. |7J the decor- 
relation bandwidth is 6Ad ec /A ~ 1%. Since the K s passband is 



6 For a given physical screen characterised by the column density 
Nl(xi,yi), Rdiff varies with A b,s , as shown in Eq. (fJJ. 

7 Nevertheless, one should consider that the gaseous screen is a 
non-dispersive medium for optical wavelengths (dielectric medium with 
an index independent of A), unlike the radioastronomy case (plasma), 
which may result in weaker chromaticity sensitivity. This specificity 



6A/A ~ 0.1 ~ 10 x SAd ec /A, a first-order estimation of the modu- 
lation index for a K s passband pattern is the value obtained when 
adding ten decorrelated patterns of speckle i.e. one multiplied 
by ~ Vl/10. The modulation index found with our simulation 
(~ 55%) has the correct order of magnitude. 

By contrast, structures of size of R re j are much less sensitive 
to the variations in A, according to the combination of ([3]i and 
©: 



ARref(A) 1 AA 1 

— = - x — = - x 10% ~ 2%. 

R ref {A) 5 A 5 



(19) 



As a consequence of the large-scale smearing of the point-source 
pattern when considering an extended source, there is no sig- 
nificant difference between monochromatic and polychromatic 
patterns for such an extended source, as shown in Fig. [7] (right). 
Therefore, in the following we ignore the impact of the Rdiff 
structures. 



3.3. Simulation of light curves 

What we observe with a single telescope is not the 2D illumi- 
nation pattern but a light curve. Because of the relative motions, 
the telescope sweeps a ID section of the pattern, at a constant 
velocity as long as parallax can be neglected. We therefore sim- 
ulated light versus time curves by sampling the 2D pattern pixels 
along straight lines, with the relative speed of the telescope. 



3.4. Computing limitations 

We adopted the same number of pixels NxN and the same pixel 
scale Ai to numerically describe both the screen's and the ob- 
server's planes; indeed, the light emerging from the screen is es- 
sentially contained in its shadow, and this choice optimises the 
screen filling. Since the two planes are conjugated, the following 
relation arises: 



NA\ = 2nR z F . 



(20) 



Because of this relation, N has to be as large as possible to min- 
imize Ai and also to get a wide useful area. This area is indeed 
restricted by the definition of the fiducial domain, which can be 
heavily reduced when simulating the illumination from an ex- 
tended source. Also a large fiducial domain is essential when 
simulating long-duration light curves. As a consequence, mem- 
ory limitations affect the maximum size of the screen and illu- 
mination 2D patterns. In the present paper, we were limited to 
patterns of 20, 000 x 20, 000 pixels. 

4. Observables 

The observable parameters of the scintillation process are 
the modulation index and the modulation's characteristic time 
scales. The main observable we used in our data analysis 
(IHabibi et al.l 201 lb) is the modulation index. It is defined as 
the flux dispersion 07 divided by the mean flux I: m - crj/I. 
We first compare the effective modulation index of simulated 
screens with the theoretical expectations, then examine the pre- 
cision we can reach on m from a light curve, and give a numerical 
example. We also show how we can connect the observed mod- 
ulation index to the geometrical parameters (Rdiff(A), A, R s etc.) 
through simulation. 



deserves more detailed studies, which are beyond the scope of this pa- 
per, considering the minor impact of the bandwidth compared to the 
smearing produced by the source size. 
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Figure 7. Simulated illumination maps (20, 000 X 20, 000 pixel of '22.6 km side) produced on Earth by a source located at zo + Z\ 
1,18 kpc through a refracting cloud assumed to be at zo — 160 pc with a turbulence parameter Rdiff (2.16pm) = 100 km. Here 
R ref (2.\6pm) = 100 000 fern. 

Top-left and middle: illumination produced at A — 2.16pm from a point-source with a zoomed detail; the contrast is 100%. The grey 
scale ranges from to A times the mean intensity. 

Top-right: the same from a K0V star (r s — 0.85R Q , My — 5.9, at \A8kpc V — 16.3). The circle shows the projection of the stellar 
disk (Rs — r s X Zo/zi). Here the modulation index is only 3.3%, and the grey scale ranges from +20% around the mean intensity. 
The bottom maps are the illuminations in K s wide band (A centra i — 2.\62pm, A/1 = 0.215pm), using the same grey scales as above. 
The modulation index is 55% for the point-source (left and centre) and 3.3% for the extended source ( right). 

The two parallel straight lines show the sections sampled by two observers located about 10000 km apart, when the screen moves 
with the transverse velocity Vj. 



4.1. Modulation Index 

For a point-like source in the strong scintillation regime the mod- 
ulation index m » 1 . For an extended source (radius r s , projected 
radius R s = r s .zo/zi) in the same regime, we always have m < 1, 
and |Narayan| ( 1 992) showed that when the small-scale (diffrac- 
tive) speckle is completely smeared, 



(21) 



in the case of the Kolmogorov turbulence. Here, 6 re f = R re f/zo 
is the angular refraction radius and 6 S the source angular ra- 
dius. In this expression, Narayan| (1992) assumed that zo « Zi, 
therefore zo+Zi ~ Z\. We use the following expression, which is 
formally identical to the previous one for zo « Z\ but can also 
be used when zo is not negligible: 



R diff 


1/3 


e ref 


7/6 


Rf _ 




2n6 s 





R diff 


1/3 


R ref 


7/6 


Rf _ 




2kR s 





(22) 



Narayan uses i 



which equals d re f/2n. 



0.035 



Ipm 



lkpc 



R, 



diff 



1000km 



r s /zi 



R Q /10kpc 



(23) 



This relation can be qualitatively justified by noticing firstly 
that Rdiff quantity has to be considered relatively to Rf when 
considering its impact on diffraction, and secondly that the con- 
volution of the point-source pattern (see Fig. |7]-up) with the 
projected stellar profile (following expression (TTTt ) makes the 
contrast decrease when the number of ^ rc /-size domains within 
Rs - rszolzi increases. The "exotic" exponents are related 
to the Kolmogorov turbulence. It should be emphasised that 
this relation assumes the scintillation is strong and quenched 
(R s /R F >R F /R diff >l). 

We checked this relation by performing series of simulations 
with different phase screens and stellar radii. We generated series 
of screens with Rdiff = 50 km to 500 km by steps of 50 km. For 
each screen, we considered different sources -at the same geo- 
metrical distances- with radii from 0.25R e to 1 .5R by steps of 
0.25R Q and computed corresponding illumination pattern reali- 
sations. The modulation indices were estimated within the fidu- 
cial zone for each 2D illumination pattern. In figure [8] the thus 
estimated m for each generated pattern are plotted as a function 
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4.2. Information from the light curves 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 

Theoretical Estimate 



Figure 8. The effective intensity modulation index m — o-j/I for 
simulated scintillating stellar illumination patterns as a func- 
tion of the theoretical modulation index. The blue dots show 
the effective modulation index values for different realisations of 
the phase screens. The red squares represent the mean value of 
the effective modulation indices with the same theoretical value. 
This plot shows the expected agreement for m < 0.15, a domain 
where the simulated screens are large enough not to suffer from 
statistical biases ( see text). 



of the theoretical value expected from expression (123V We note 
that the modulation has relatively large scatter. This is because 
the fiducial zone is not large enough, and it contains a limited 
number of regions with sizes of R re f (the scale that dominates 
the light variations). In some cases there are very few distinct 
dark/luminous regions (as in figure [5}. The number of such re- 
gions within the fiducial domain is Nr ~ d 2 /R* e p where d is the 
size of the fiducial zone. After discarding the cases with Nr < 
5, we computed the mean and root mean square of the effective 
m values, for each series of configurations with the same the- 
oretical modulation index. Figure [8] shows that relation (j23j is 
satisfied for a modulation index smaller than 0.15. The higher 
values of m are systematically underestimated in our simulation 
owing to the limitation of the screen size, which has a stronger 
impact on the number of large dark/luminous regions in this part 
of the figure. This number is indeed smaller (though it is larger 
than 5), and the chances for sampling the deepest valleys and the 
highest peaks are fewer, therefore biasing the modulation index 
towards low values. 

Therefore, by measuring the modulation index from an ob- 
served light curve and using an estimate of star type (i.e. radius) 
and distance zo+Z\ ~ Zi , we can constrain Rdiff(A) from this re- 
lation d23lFl even with poor knowledge of the screen's distance 
Zo, considering the slow dependence with this parameter. This 
technique allowed us to infer constraints in Habibi et al. (201 lb) 
on the gas turbulence within galactic nebulae (zo ~ 80 - 190 pc 
and Zi = 8kpc), and upper limits on hidden turbulent gas within 
the galactic halo (assuming zo ~ lOkpc and zi + Zo = 62kpc). 




Figure 9. Light curves extracted along the 3 horizontal white 
lines for two illumination patterns. Left column: 2D pattern from 
a point-like source in K s band with Rsff = 300 km and R re f « 
28000 km. The modulation indices of the three light curves differ 
by less than 5% from the 2D pattern modulation index. Right col- 
umn: the illumination pattern for an extended source with R s * 
41000 km, through the same refractive screen. The modulation 
indices fluctuate by more than 30% around the 2D pattern in- 
dex, implying the necessity of longer light curves for a better 
statistical representativity. The distance scale is common to both 
patterns. The circle shows the projected star disk. The 3 light 
curves from the right column are not completely decorrelated, 
because of their common proximity to the same large positive 
fluctuation. 



If a light curve is long enough (or equivalently if the obser- 
vation time is long enough), its series of light measurements rep- 
resents an unbiased subsample of the 2D pattern. For instance, 
the left hand panels of figure |9]represent the illumination pattern 
of a point-like source and three associated light curves. Here, 



R 



ref 



28, 000 km and the modulation index m 



poml 



= 1.18. The 



9 In dHabibi et al. 1 201 lb), we used a simplified relation, with no sig- 
nificant impact on the resulting constraints. 



light curves are extracted from three horizontal parallel lines 
with lengths of ~ 3.5 x 10 5 km. The corresponding time scale de- 
pends on the relative transverse velocity. The lines are selected 
far from each other in order not to be affected by the fluctu- 
ations of the same regions. Modulation indices along the light 
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r(km) 



r(km) 



Figure 10. Simulation of illumination scintillation patterns with associ- 
ated light curves from an extended source assuming a refractive screen 
with non-Kolmogorov turbulences (left B — 3.1, right B — 3.9, see text). 
Here Rdiff = 100 km, R re f «8x 10 4 km and the projected star radius is 
Rs ~ 036R ref = 28,800/tm. 



3 
o 

Q. 



|3 = 3.9 
P = 3.1 



f)- 1 I 

"diff ^ 



log (q (1/km) ) 




0.06 0.08 0.1 0.12 0.14 

Kolmogorov Modulation Index 

Figure 11. Top: Two different phase screen power 
laws. Bottom: The corresponding modulation indices 
as a function of the expected modulation index for 
the Kolmogorov turbulence. The 3 indices plotted at 
a given abscissa correspond to screens with B — 3.1 
(red), 3.67 (Kolmogorov, blue, along the diagonal), 
and 3.9 (black) with the same Rdiff- As the power 
law gets steeper, a larger modulation is expected. Blue 
dots show the Kolmogorov turbulence case. 



curves differ from the 2D's by less than 5%. When the mod- 
ulation is characterised by both R re f and Rdiff, a light curve 
that spans a few R re f is a sample of the 2D screen, which is 
large enough to provide a good approximation of the scintilla- 
tion modulation index for a point-like source. The right hand 
panels of figure [9] show the 2D pattern for an extended source 
with R s « 41,000A77i. Here, R s > R re f and the flux fluctuations 
are smoothed, characterised by the unique length scale R re f, 
and have a much smaller modulation index m exten d e d = 0.04. 
The light curves are extracted in the same way as the point-like 
source within the corresponding restricted fiducial zone (see fig- 
ure 0, therefore they span ~ 2.5 x 10 s km and statistically in- 
clude a little bit less than 10 Z? re /-scale variations. Because of 
this statistically short length, the light-curve-to-light-curve es- 
timates of m exte „ded typically fluctuate by ~ 1/VlO ~ 30%. 
The fluctuations on m exten d e d estimates can only be reduced with 
longer light curves. To get a precision value of 5% on m exten d e d, 
the light curve should indeed be as long as ~ 400 x R re f- As 
an example, if we observe through a turbulent gaseous core 
with Rdiff = 200 km in B68 nebula located at zo = 80 pc at 
A = 2.16 pm (R re f ~ 27, OOOfcra), and assuming VY ~ 20 km/s, 



an observing time of ~ 150 hours is needed to measure the mod- 
ulation index with this precision of 5%. When searching for un- 
seen turbulent media located at unknown distances from us, the 
diffraction and refraction radii are unknown, and we can only 
obtain a probability distribution of the observation time for a re- 
quested precision on the modulation index. 

The other important information carried by the light curves 
are the characteristic time scale between peaks t re f = Rref/Vr, 
associated to the refraction radius (Fig. |7]upper left), a correla- 
tion duration t s = R s /Vt, associated to the source radius (Fig. 
[7] upper right) and possibly tat/f = Rdiff /Vj, associated to the 
small diffractive speckle structure (Fig. |7]lower left); the latter 
could be detected in exceptionally favourable cases (source with 
very small angular size, assuming the strong diffractive scintil- 
lation regime) with a powerful detection setup such as LSSlFl 



10 The detection condition would be rjzi % 10 x Rdiff /zo (f° r 
which the projected stellar disk includes less than ~ 100 speckle spots), 
assuming a setup able to sample the target star with < 1% photo- 
metric precision every few seconds. The scintillation of an AO type 
star (R = 2AR Q ) in the LMC (magnitude V = 19.4) seen through a 
screen with Rdiff = lOOfon at distance 30 pc is an example of such 
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The time power spectrum of the expected stochastic light curve 
should show a slope break around frequency and possi- 



bly around t d L. ( |Goodman & Narayan 1985) , which should al 



low one to distinguish it from purely randomly fluctuating light 
curves due to photometric noise, and to extract constraints on the 
scintillation configuration. Estimating t s would be challenging, 
since the extended source acts as a low passband filter on the 
point-source pattern; therefore, the imprint of the projected ra- 
dius is essentially within the attenuation of the light curve auto- 



correlation (Goodman & Narayan 2006), and more specifically 
-as mentioned above- within the modulation index, that can be 
polluted by many observational artefacts. These potentialities of 
the time studies need more investigation and should be discussed 
in more detail in a forthcoming paper. 



5. Probing various turbulence laws 

Up to now, we have focussed on the standard Kolmogorov tur- 
bulence. In this section, we investigate a possible deviation from 
the Kolmogorov turbulence lavvFl We chose two different phase 
spectra with f3 = 3.1 and f3 = 3.9 in relation (|4j. To study the 
corresponding scintillation modes, we generated two series of 
phase screens according to the spectra and computed the illu- 
mination patterns from an extended source with R s ~ 0.36 R re j 
through each of them. The patterns are represented in figure [TO] 
with corresponding light curve samples. For both images, Rdiff 
= 100 km, R F = 1150 km, and R ref ~ 8 x 10 4 km. The pat- 
tern with p = 3. 1 (left) shows a small modulation index m(3.1) = 
0.04, while the pattern with /3 = 3.9 (right) shows a much larger 
modulation index ra(3.9) = 0.22. As can be seen from its visual 
aspect, the turbulence with a higher exponent produces stronger 
contrast on large scales compared to the other one (/? = 3.1). To 
understand the origin of the difference, we compared the two 
phase spectra at the top of figure [TT] The steeper spectrum (J3 = 
3.9) has more power for fluctuations on large scales. Moreover, 
by integrating equation (0), we computed that the total power 
distributed from R re j to Rdiff is about an order of magnitude 
larger for /3 = 3.9 than for /3 = 3. 1 . Therefore, the larger the /3 ex- 
ponent is the stronger flux fluctuations are produced. As a con- 
clusion, the detection of scintillation should be easier for turbu- 
lences with steeper spectrum, because a larger modulation index 
is expected. This is illustrated at the bottom of figure QT| where 
we plot the modulation indices produced by three screens with 
different fj values (including the Kolmogorov case), as a function 
of the "classical" Kolmogorov turbulence expected index. In this 
plot, each abscissa corresponds to a distinct Rdiff value, which 
is common to the three screens. By increasing /?, we increase the 
modulation, especially for high m values. 



6. Discussion: guidelines provided by the 
simulation 



The observations analysed in (IHabibi et al.l 201 lb) were in- 
trepreted by using our simulation pipeline. But we have used 
our simulation not only to establish connections between the ob- 
served light curves and the scintillation configuration, but also to 
define observing strategies as follows. 

Firstly, a correct sensitivity to the scintillation needs the ability 



a favourable configuration, which could be discovered by the LSST 
JLSST Science Book 20091 . 

11 For the observed supersonic turbulence (with /? > 11/3) see 
|Larson (1981)1 



to sample, with < 1 % photometric precision at a sub-minute rate, 
the light curves of small distant background stars (M ~ 20-21), 
which have a projected radius small enough to allow for a mod- 
ulation index of a few percent (typically < R at lOkpc). Our 
study of the time coherency shows that the usual large passband 
filters can be used without significant loss of modulation index. 
Since the optical depth of the process is unknown, a large field 
of view seems necessary for the exploratory observations, ei- 
ther toward extragalactic stellar sources within LMC or SMC or 
through known gaseous nebulae. To summarise, an ideal setup 
for searching for scintillation with series of sub-minute expo- 
sures would be a ~ 4 m class telescope equipped with a fast read- 
out wide field camera and a standard filter (optical passband to 
search for invisible gas towards extragalactic sources, infrared to 
observe stars through visible dusty nebulae). 
Secondly, our work on the simulation provides us with a guide- 
line to find an undisputable signature of scintillation. The first 
possibility consists in the search for chromatic effects. Subtle 
chromatic effects betwen the different regimes (associated to 
different time scales) have been shown, but will probably be 
hard to observe. Figure [7J (left and centre) shows the expected 
speckle image from a point source. The position/size of the 
small speckles (characterised by RdiffW) are very sensitive to 
the wavelength, and it is clear that a desynchronisation of the 
maxima would be expected when observing such an idealised 
point source through different (narrow) passbands. But that real 
sources have a much larger projected radius than the speckle 
size completely screens this chromatic effect. The impact of the 
wavelength on the position/size of the wide (refractive) spots 
(Fig- Gkight) is much weaker (following X ' according to the 
combination of expressions © and (|9]i)- Therefore, only a weak 
chromatic effect is expected from an extended source, even when 
observing with two very different passbands. For a clear sig- 
nature of scintillation, it seems easier to take advantage of the 
rapidly varying luminosity with the observer's position within 
the illumination pattern. As a consequence, simultaneous obser- 
vations with two ~ 4 m class telescopes at a large separation (few 
10 3 km) would sample different regions of the illumination pat- 
terns (see Fig. [7J bottom right), and therefore measure (at least 
partially) decorrelated light curves as shown in Fig. |9] (right). 
The decorrelation will be complete if the distance between the 
two telescopes is greater than 2 x max(R re f, R s ). A single obser- 
vation of such a decorrelation will be sufficient to definitely con- 
firm the discovery of a propagation effect that cannot be mim- 
icked by an intrinsic variability. 

7. Conclusions and perspectives 

Through this work, we have simulated the phase delay induced 
by a turbulent refractive medium on the propagation of a wave 
front. We discussed the computational limitations of sampling 
the phase spectrum and to obtaining large enough illumination 
patterns. These limitations will be overriden in the near future 
with increasing computing capabilities. The illumination pattern 
on the observer's plane has been computed for the promising 
strong regime of scintillation, and the effects of the source spa- 
tial and time coherencies have been included. We have estab- 
lished the connection between the modulation index (as an ob- 
servable) of the illumination pattern with the geometrical param- 
eters of the source and the strength of the turbulence (quantified 
by Rdiff)- Furthermore, we showed that when the spectral index 
of the turbulence increases (as in the case of supersonic turbu- 
lence), the detection of scintillating light curves should be easier. 
These simulation studies and, more specifically, the modula- 
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tion index topic were successfully used in our companion paper 
dHabibi et al.1 201 lb) to interpret our light curve test observations 
of stars located behind known galactic nebulae and of stars from 
the Small Magellanic Cloud, in the search for hypothetical cold 
molecular halo clouds. 

Time scales, such as t re j = R r ef/Vr and t s = R s /Vt (where VY is 
the relative velocity between the cloud and the line of sight), are 
observables that we plan to study further and in detail. Their ex- 
traction can be done through analysing the time power spectrum 
of the light curves, and it should give valuable information on the 
geometrical configuration, as well as on the turbulent medium. 

The observing strategy has been refined through the use of 
the simulation, and we showed that observing desynchronised 
light curves simultaneously measured by two distant four-metre 
class telescopes would provide an unambiguous signature of 
scintillation as a propagation effect. 
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